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PROGRAM FOR THE DESIGN OF AN AXIAL COMPRESSOR STAGE 
BASED ON THE RADIAL EQUILIBRIUM EQUATIONS 


by 
Kyriacos D. Papailiou 


i. introduction 


The present report describes a program which has been established 
for the design of an axial flow compressor intermediate stage using the 
radial equilibrium equations as described in ref. 1 (pages 439 - 45h). 
An intermediate stage is defined as one for which inlet and outlet 
velocities are essentially the same. 

The flow is assumed compressible and axisymmetric. The inlet 
conditions, some of which are data to the program, are assumed to have 
been produced in preceding stages and are referred to ambient condi- 
tions). Energy gradients at the inlet, which may be caused by 
non-uniform energy additionin previous stages,as well as entropy gradients 
at the inlet, which may be due to non-uniform energy dissipation in 
previous stages, are taken into account. 

The losses and the corresponding entropy increase through the rotor 
and the stator are calculated according to reference 3, and non-uniform 
energy increase in the rotor can be specified. 

The inner and outer walls may have arbitrary shape. However, when 
the curvature of the streamlines (in the way described in reference 1.) 
is taken into account, then the inner and outer walls, although maybe 
tapered, have to be straight. 

The problem is considered from the designer's point of view. Con- 
sequently, specified will be the work distribution, and the quantity 
Vuz+Vue 


Uj tu 
The quantity 2 gives the theoretical reaction factor for an inter- 
at 


will be specified along the radius (see symbol table for symbols). 


mediate stage with constant axial velocity for incompressible flow. This 


quantity was chosen instead of the actual reaction factor as the expression 


fia 7 THERMODYNAMIC DIAGRAM 
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of the actual reaction factor is rather involved in the general case 
(see Appendix A). The thermodynamic diagram of the process through 
the stage is given in Figure (1) while the general layout is given 
in Figure (2). 

In the following the theory will first be developed and then 
the program will be described. Symbols are defined in the symbol 
table and the FORTRAN Symbol Table. 


SR i 
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FIG. 2. GENERAL LAYOUT 


2. General Equations Used (See Ref. 1) 


Vat ry ace (1) 
dR a 
where F = se ats = = a S (2) 
: a ee eta as (3) 
with ace a = — == 2K (4) 


The minus sign holds for stations (e) and (2) and the plus 
sign for station (1)* 


¥Note that there is a discrepancy between the way equation (4) is inter- 
preted here and in ref. 1. However, positive curvature is considered here 
in accordance with the way equations (1), (2) and (3) have been developed. 
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pS = =tAR 
where tan A, = tan \, = tan No = oa (5) 
L,+ L 
i= ao (See Fig. 2) (6) 


K takes the values 4 to 6. 


The solution of equation (1) is given as follows 


2 z Bae R’ 
V, exp (. \) ran) Ya, - | cox (+ | rar’) dR 
Rn Rn Rh 


R R R! 
_fFdR _/FAR 3 f Far’ 
= ve 2 hae e Rh . " Ge ®h dR (7) 


3. Dimensional and Non-Dimensional Quantities 


The basic quantities used are expressed as follows 
Lengths in (ft) 
velocities in (ft/s) 
angular velocity w in (rad/s) 
ft- fteib) 
slug 
ESN AS ye) 
slug,°R ane 3°R 


enthalpies in (£t2/s°) or ( 


entropies in ( 


the specific heat Cp in ——- or (—— = 2a) 


the density in (slug/ ft?) 


Consequently F has dimensions (ft71) 


G has dimensions (£t/s°) 


Calculations will be formed with non-dimensional quantities. As 
reference quantities will be used: The angular speed ”, the mean radius 
Rn» where 


R 
$1 + Rhy 
i= on ae (7a) 


the atmospheric pressure, temperature and density. We shall call all 
non-dimensional quantities starred quantities and denote them with a 
star. 


- Lengths will be non-dimensional over R,. Consequently 


Aas 
Rn 
i 
29) =} 
Loa 


- Velocities will be non-dimensionalized over WRy,. Consequently 
(where Vrer = Ry) 


Enthalpies will be non-dimensionalized over uPR * 


h H 


h*¥ =———_ 5>—- HE = ———— 
oP Rye Ry 


Entropies will be non-dimensionalized over CS 


Densities will be non-dimensionalized over the atmospheric density 


- Pressures will be non-dimensionalized over the atmospheric pressure 


4, Problem Formulation 


Using non-dimensional quantities (in the described way), we arrive 


at the following equations 


a(V_* 
—sar «+ * (V,*) + GF = 0 (1a) 
where 
OW,” 30(S/Cp) 
pee Perea Tesppeeues O64 Gy) (2a) 


Cx = 2 Yu* o(R*V)*) a 2oH* + IC 2nRe y= (v,*)) a(S/Cp) (3a) 


Re OR OR* OR* 
and 
Fx = F Rm 
(8) 
G* = G/weR, 


We shall consider now Figure 3. Considering stations @ (2) and (3) 


FiguRE 3. 


we can see that the angle \ can be found from 


bOARS HSR3I=F RI 
tanh = oF BL 


and still the curvature be calculated using §R at station (2) which is 
now negative measured from the straight line @) - ey instead from the 


straight line (e) - ey The curvature expression remains the same if 
this sign change is taken into account. 


5. Necessary Inlet Conditions 


Quantities will be specified at equal radial distances at the inlet 
The distribution of total enthalpy H and the distribution of entropy s 
must be specified at the inlet. This distribution of entropy am total 
enthalpy will be considered to have been developed through an adiabatic 
process (previously existing stages), starting from uniform atmospheric 
conditions. The atmospheric total pressure and temperature must be also 
given. Entropies will be measured considering as origin the atmospheric 
entropy. 
Considering an adiabatic compression 
from (A) (atmospheric conditions) to 
(@) (conditions at the stage inlet) 
along a stream surface, we have 


(see Fig. 4 for symbols) 


HR ” 
Sh - Spx = Ce In Soe = Sp- S,=AS 


Hp* = Hp @48/ Cy 


F/G. 4. Ad1IABATIC COoMPRESS/ON 
W/TH FRICTION. 


from where Hp* can be calculated. For an isentropic process we have 
also (from state A to state B*) 


-As/C C= 
Pep* Pep HpY W/(v-1) [Hee p \(y-1 
Pea Pea Ha Cy Ten 


from where the total pressure at the compressor inlet can be calculated. 


The stagnation density Py. is given then as 


geste. kept 
‘= = 


6. General Method of Solution 


As can be seen an iterative process is necessary to solve the 
problem. A first approximation of the streamline position is assumed 
and on this basis the radial distribution of the axial velocity vo 
is calculated at stations (1), (2) and (3), that satisfies the 
continuity equation on the whole. 

Then, a new approximation of the streamline position is achieved 
by requiring that the same mass flow passes through individual stream- 
tubes at stations (2) and (3), the radial positions at station (1) 
remaining the same. 

A new V,*-distribution is then calculated and the iteration continues 
until the specified number of iterations is achieved. The printed 
error of each iteration gives us an indication of the convergence of the 
procedure. 

Having established the equations to be solved, a general flow diagram 


is given in Table I describing the general layout of the program. 


7. Detailed Calculations (Following flow diagram of Table I). 


Block G) 
At stations G) () ana(3) - 


equidistant radii are considered 


FIG. 5. 


Block G) 


The total enthalpy increase inside the rotor is 
AH = 8 (Ro Vup - Ry Vuy) 


* = = % Te es * * 
on Aah eS . Ro* Vyp* - R,* Vy 


2 a 


The reaction factor is 


Yu, + Vuo : Vy /R, + Vup/ OR, i Wi Vor 
UL + Us, U,/oR,, + Up / WRy Ry* + R,* 


RF = 


or finally 
* a . %& * 
she RF(R,* + R,*) - AH*/R, 6 
u - & oa 
a Des [R., 
RF(R,* as R,*) + AH*/ Be 
U 
e aa ae 


Vu3* = Vay" 


Block G) 


Calculations are performed without losses in the stage nor curvature 
effects for the first iteration. Corresponding terms are set here 


equal to zero. 


Block( 4) 


The curvature term is calculated according to equation (4) rewritten 


here 


Se ee tio a (11) 


Note that in the way Vavra describes the curvature in page 453 of ref. (1), 
there exists an inconsistency in sign with the derived equations of motion. 
In fact the plus sign belongs to station z = O of ref. (1) or stations (1) 
and (3) here, while the minus sign to station (2) here (if 6R is taken to 


be a positive quantity). The expression for the curvature term then becomes 


2 OVy 5R 
CRTERM =- — Ree ek —— ie 
Va Zn — Te 
x * 
' Scene td 
Ke 
2 ge tee ie ae tet ~ E (12) 

= Lg" +Lp Hi = & +LR ‘i 

. a lose Tae 


In this expression 6R* is positive in the way shown in fig. 3 and (+) sign 


applies to stations (1) and (3) while (-) applies to station (2). 


The losses in total pressure are calculated in a subroutine, the 
theoretical basis of which is given in appendix B. This was decided in 
order to have the freedom of introducing any loss-correlation we desire. 

The subroutine will furnish to us the decrease in total pressure non- 
dimensionalized over the atmospheric pressure. 


The resulting entropy increase will then be calculated as 


where P, corresponds to the total pressure level without losses. For the 
rotor this total pressure is Pea (see fig. 1) which is the corresponding 
total pressure that would result if the addition of enthalpy Ho- H, were 


al 
done isentropically. Then we have 


! ¥ noe * i 
St aif Heh yen aa Uae ale, Nya) (43) 
Dee Hy ae H,* 
and 
y 
* 
a = veal 
= = * 
85-8, __ Rg _ 1" \Hy* (AP,*) 3 : 
C Cy Ho” ae 
* - 
rea Hy] 
; o Na 
APE [Hoe \ (YoU) 
ees Sa lt = i = Sak Seger ( 
=- -=— | => = - 14) 
Cp Pee ee 2 a 


For the stator IE comes to be equivalent to Pio and the increase in 


entropy from station (2) to station (3) is given as 


Pe. (Ae) 
@--#n (= = 


C Se Pte 
or 
as (AP* 
3352 Rg ( = 
Ein Pett, Tapes (a ute eet In {1 - ———jJ = S,* - S,* (15) 


2UG) 


Block G) 


The subroutine DERIR (see description in ref. (4)) is used to find 


the derivatives of a function given at discrete points. Thus a ; — 
3(R¥Vy* ) : : 
—spe:~C~«CS': calculated. Then the functions F* and G* are calcu- 


lated as given in equations (2a) and (3a). 


Block(6) 


The integrals of equation (7) are calculated (FUNC1 and FUNC2) and 


the distribution of (ies is considered in the following form 
(Vv se = FUNCL (Vv xe + FUNC2 
a ah 


An iteration is initiated using as starting Vay,* Value one which does not 


make the quantity Gee negative. The Seer scheme is described in 
ref. 4 (iteration subroutine ITERN). 

For each value of (Vap*)> the es) - distribution is calculated 
for the station under consideration. The iteration process ends when the 


continuity equation is satisfied. 


Block ©) 


The continuity equation is considered in the following way. Consider 
the actual axial velocity distribution at a section, V,* = Va*(R*). Then 


the corresponding mass flow is 


oR 


aan 
real Rh 


Q* OTIPxV .* R¥AR* (16) 
Consider now a hypothetical situation where the part of the inviscid flow 
distribution of ee is extended inside the wall boundary layer regions 
up to the walls. Then using this distribution we get a mass flow rate Q*, 


where 


Ry % 
Qe = [  ” anpxv R¥aR* 
Rn ainviscid 
then 
> ae 
@ 8 real 
x *) qi f * insi 
and Ve" viscid (R*) differ from V,*(R*) only inside the wall boundary layer 
regions. 


ol 


We can say now that normally a flow rate Q* would pass through the 
area given,if it were not for the wall boundary layer presence. Defining 


a blockage factor Ky as 


Q* real 
Ks Eee (Gia) 
Qx 
then 
Q* Ret 
“real _ Qe = | Onp*V_* R* aR* (18) 
R ae BO ie 
K. h inviscid 


where the mass flow rate has been non-dimensionalized as follows 


Q 


Patm 


Qe = 
Rn? 
In the program the inviscid velocity distribution Vg* . ., is calculated 
inviscid 
and throughout the whole report this is being referred to as Vat The 
mass flow calculated at each station then is Q* and then Ky, Q* is compared 
5 % s 

WEEE real 

The calculation of the density p* needed in the continuity equation 
is performed as follows: The absolute velocity is calculated as 


poo 
) je 


(Vy*)® + (vy*)* + (Vy 


U(ve)® + (v,*)® + (v,*)® tanPa (19) 


The static temperature is calculated as 


Ve 


He (vx)® 2 2 


eee ee (20) 


The velocity of sound is calculated as 


ge gs, ERTL (21) 


Nive (22) 
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Then the density is: 


al 
p Co zy Wig) Lory sAdat 
ease aes tat bP a 
atm “g ~ Patm 
i 
* yor: 
a Cy Pe P atm ( us Y-1 #) y ) 
wreRe 2 2 
& H* w Ry Patm 
aD 
P % -(Y-1) 
t Y-1 we 
= — 4H* 1+— 2 
wie 5 (23) 


The calculation of the mass flow rate at each radial position 
is performed with the subroutine SUMAT which uses Simpson's rule for 
unequal intervals. The subroutine SUMAT makes use of the subroutine 
INTPO) for interpolation which is described in detail in ref. 2. 

Once these calculations are performed the program gives control 
to the subroutine ITERN which performs the iteration, until convergence 


has been realized and the continuity equation is satisfied. 


Block (8) 


The radial position of the streamlines at station (1) is considered 
the same throughout the calculation and in this step the radial positions 
of the streamlines at stations (2) and (3) are found with the condition 
that the same mass flow is allowed to pass through each streamtube with 


the newly calculated axial velocity distribution. 


Block (9) 


The square root of the sum of the squares of the difference between 
the old and new position of the streamlines is calculated for stations (2) 
and (3). This is considered here as an indication of the convergence of 
the procedure. Additionally in this step all pertinent quantities that 
have not been calculated up to now are calculated (static and total pressures, 


relative velocities and mach numbers, angles, etc). 


13 


Block Gg 


Once the prescribed number of iterations has been performed, 
some additional pertinent dimensional and non-dimensional quantities 
are calculated (a description of the calculation is given in Appendix C), 
and all the results are printed. 

The whole procedure is executed as many times as additional set 


of data exist. 


8. Description of the Use of the Program 


This program has been constructed as a complement to ref. 2. An 
effort was made to take into account entropy and energy gradients (entropy 
gradients existing at the inlet or introduced through the loss correlations 
after each row,and mergy gradients existing at the inlet or introduced by 
a non-uniform work distribution in the rotor) and compressibility effects. 

The non-dimensionsalization proposed in ref. 2 which does not take 
these effects into account is rendered thus incomplete and it was decided 
to carry out the calculation in the non-dimensional form proposed in 
reference 2 and modified slightly as described already, introducing, however, 
the data in dimensional form, reflecting thus the Mach number level. 

A detailed description of the program has already been given. The 
meaning of the weight factor has to be explained here. It happens sometimes 
in complicated cases that the iteration procedure diverges when the calculated 
corrections for the new streamline position are used in the whole. If, however, 
a fraction of the corrections is considered and introduced for the next itera- 
tion,the iteration procedure may be forced to converge. The weight factor 
( 0. < W.F. < 1.) introduced as data to the program defines the fraction of 
the correction to be used for the following iteration loop. 

Cases have been already run for no curvature effects and no losses. 

Then losses were introduced and curvature effects. It was found that to 
take into account the curvature effects and have an converging iteration 


process, a value of the weight factors smaller than unity ought to be used. 
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To facilitate the use of the program a table with a typical input 
has been prepared and given in Table II. The maximum number of stream- 
lines is taken to be eleven. The listing of the program along with the 
subroutines in use is given in Table III. The results of the already given 
typical input are given in Table IV. The explanation of the symbols 
used along with all the pertinent parameters used in the program are given 
in the FORTRAN Symbol Table. The output symbols not described in the 
FORTRAN Symbol Table will be found in the Output Symbol Table. For 
dimensional quantities the already given dimensions are used. For non- 
dimensional quantities the already described non-dimensionalization has 
been done. In the FORTRAN Symbol Table the dimensional variables are 
given. The non-dimensional ones are denoted in the program with the letter 


S at the end of the name of the dimensional quantity unless otherwise stated. 


1b 
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APPENDIX A 


THE REACTION FACTOR WITH AXIAL VELOCITY VARIATIONS 


Ry # By 
AH = 0(Ro Vyo - Ry Vu; ) 


yh LPS 
APom 


Assuming no losses, for 
incompressible flow 


Ae = Fara) = 4? to) eer “Fey per, 
= -u4° ne 
welt 2 D 
= $0(W,2-Wo“) + 4 (Up°-U,") 
2 
£16. 1(@) VELOCITY TRIANGLES = 3°( Wy) - Wyo? + Vay -Vao®) 
FIGURE 1A + 30(U>*-U,*) 
AP on = W(Ro Vuo - Rj Vy) P 


which turns out to give a rather complicated formula for the reaction 
factor. Consequently, we shall assume that the quantity 


Vy + Vy. 


1 2 
Wy + Uo 


is specified along corresponding radii. 
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APPENDIX B 


The correlation described in reference 3 will be used for the 
calculation of losses (see also reference }+) 


The loss coefficient ¢ . is defined for rotors as 


ere 
Co = (Azo ee (B1) 
= = * 
Ri aes Pey* - Py 2 


where AP, is the mass averaged total pressure loss, and for stators as 


C (AP,)s (APo*)s 


P CC —_—_—_————— ———_—_———— 
8 iPus\- Po Pip” - Bo* 


(B2) 


This definition holds for incompressible flow where the density is assumed 
constant 


The diffusion factor D is given for rotors as 


Wo Ry Wu, - Ro Wao 


D-=1-—+ (B3) 
e Wy OR (Ry + Re) Wy 
and for stators as 
eo ee (Bh) 
* a: es * * 
or D=1 ie ip arias coh (B3a) 
% i ¥ = 
R W, AGC Ty Cts 
Vax Ro* Vup* - R2* Vya* 
Dp - See ne (B3b) 
2 Sayae 3 2 
then 
2°R 2,91 2.02 
Cp = <8 Bp |o.008 + 0.0639 (Dp + 0.1) + 0.228D,, 
[reaple "| (ess 
Ro, - R,, Ro,* - Ro* 
where =o SS (B6) 
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For angles By 2 hee. see Fig. 1A, a correction is offered, where 


instead of Cpe the loss coefficient (Spo) cor is considered, where 


(ecg a a coe By/2 Ly + 8) (B7) 


For stators 


eo 
Spy = —— [0-008 #rentG39 (0.4 0:1) > 410,057 De 


Soe 3 
ne er (B8) 


: é fe) 
where no correction is offered for a, > 45~, 


3 


The shock losses are calculated as follows: 


The amount of supersonic turning is 


re a (6; - Bp) for the rotor (B9) 
2 eee (0 - 03) for the stator (B10) 


The peak suction surface Mach number is then obtained from 


(Mg) = 1.095 + 0.03395 Av + 1.086 (Mp) - 1.00)1-372 = (B11) 
for a rotor 
(Mg, ,) = 1.095 + 0.03395 Av + 1.086 (Mp - wOO)\n (B12) 


or Wa) Sstacor 


The shock losses are then calculated on the basis of the mean Mach number 


M + M 
R S 
M - a for ay rotor (B13) 
M, ar Ms, 
rear for a stator (B14) 


eT 


as 
, - [ove es oH ie 
(y-1)Me+2 OvMe- (y-1 


re ae Eee a -y/(y-1) 


(B15) 


where the final loss coefficient is obtained as 
C= @+ Con 


If the inlet Mach number is smaller than unity, the Mach number Moy 
is calculated by taking the inlet Mach number MR or Mp to be unity. 
The mean Mach number M then is calculated using the actual inlet Mach 
number MR or M,. If M< 1 no shock losses are assumed to exist. If M>1l1. 
then the shock losses are calculated using equation (Bl5) and M=1k. 

As pointed out earlier, this procedure is adopted only if the inlet Mach 


numbers are smaller than unity. 
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APPENDIX C 


The following additional calculations are performed in the program. 


(a) The total pressure increase is 


* ¥ ¥ 
AP eS, oe, for the rotor (el) 
R 2 al 
¥ * * 
AP, = P. eae for the stage (c2) 
Si 3 1 
(b) The static pressure increase is 
AP. =P. -P for the rot (c3) 
Roto 7 Fy or e rotor C3 
* ¥ * 
AP. =<P. =P for the stator (ch) 
S 3 2 
AP. Ee ops for the st (C5) 
aT 3 mae or e stage 5 
(c) The efficiency is calculated as follows 
Cor m, 
Tp = = total to total efficiency for the 
rotor (c6) 
Ae Saas 
Me = a a ie static to static efficiency for the 
3 e stator (C7) 
ae - *.,? 
_— is Pee 
yp) ' =——aAGT_ total to total efficiency for 
7 the stage (c8) 
(T, - T, ) 
(Th ) ee static to static efficiency for 
le aS aa 
the stage (C9) 
The calculation of the temperatures is being done as follows 
* _* 
-(8,-8, )/Cp a(n eee) 
Ty = Te = © (C10) 
Zhe 2 2 
is 


eg 


-(S,-8,)/Cp -(S,-S 
3. tT, e - T. e (Cia) 
-(S- S,)/Cp Me -(S 


ae i 
tan 3 = 


(C12) 


(d) Mass averaged quantities are calculated. Assuming that we 


want to calculate the mass averaged value of the quantity Y, then 


Reap fetes 
* 
|) 2apR eetowan ON) ene nue eae 
7 aD as 
Yer— igs Nore: > 2 (Ci3) 
Qea/X 
(one x 
tip 
en DESC ts caer ie 
=* R V, Y dR 
= i Bee (c13a) 
a ub 
-eal 
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SYMBOLS 


absolute velocity 
relative velocity 
peripheral velocity 
total enthalpy 

static enthalpy 

total pressure 

static pressure 

entropy 

radial distance from compressor axis 
angular velocity 

axial length (see Fig. 2) 


angle between axial direction and streamline direction 
at stations (1), (2), and (3). 


radial distances defined in Fig 2. 

functions defined in equations (1), (2), and (3) 
specific heat at constant pressure 

axial distance 


constant used for the evaluation of curvature effects 
(see equation (4)). 


static density 

stagnation density 

absolute angle measured from the axial direction 
relative angle measured from the axial direction 
gas constant 


ratio of the specific heats at constant pressure and constant 
volume 


3h 


Ny static temperature 


Ty, total temperature 

RF reaction factor 

AH total enthalpy increase in the rotor 
AEs loss in total pressure 

kK wall boundary layer blockage factor 
Q mass flow rate 

a velocity of sound 

M Mach number 

AP static pressure increase 

AP, total pressure increase 

D diffusion factor 

o solidity (chord to pitch ratio) 


AR,Ag non-dimensional quantities defined in equation 


1 efficiency 

Es profile loss coefficient 

sh shock loss coefficient 

Ap amount of supersonic turning 
M average mach number 

(g total loss coefficient 
Subscripts 

al station ahead of the rotor 

2 station behind the rotor and ahead of the stator 
3 station behind the stator 

a axial direction 

u peripheral direction 
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is 


ref 


cor 


rel 


radial direction 
isentropic 
entrance (see Fig. 2) 
rotor 

stator 

hub 

tip 

mean 

atmospheric 
reference 

stage 

suction surface 
corrected 


relative 
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Start 


Read required Data 
Set up a first approximation of the streamlines Gi 


From distribution of R.F.and AH* calculate V8, Vu, ,VU; 


fochisdhetire: loops Assume no losses in rotor and ©) 
9 stator and no curvature 
; effects 
Calculate Curvature terms,losses and entropy distribution 


Calculate entropy and enthalpy gradients and fromthem the 6) 
functions F°and G* 


‘ Site : : : e 
Solve radial equilibrium equation assuming a value for Va; (6) 


Test continuity @ 


Select new |NO ae — 
value of Va; Is pony Satisfied 
“w Yes 


Calculate new values of streamline positions (R32, R's) |-(8) 


alculate error for positions (2) and (3),print it and S) 
calculate additional pertinent quantities. 


has number 
of iterations pre- 
scribed inthe pro- 
gram realized 2? 


NO 


YES 


Calculate additional 
quantities and print (10) 
results 


iba 


NOTE: This procedure is repeated as many times as necessary to 
cover all cases introduced as data along with the program 
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FORTRAN IV G LEVEL 
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NON DIMENSIONALIZED QUANTITIES 
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CALCULATION OF CURVATURE TERM AND ROTOR AND STATOR LOSSES (STEP 4) 
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